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decompose the tomographic data. The wavelet packet is selected from a dictionary of wavelet packets using the best basis algorithm. 
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METHOD FOR TOMOGRAPHIC RECONSTRUCTION WITH NON- 
LINEAR DIAGONAL ESTIMATORS 

SPECIFICATION 

PPT atpd APPTJCATION 
5 This application claims priority from U.S provisional application Nos. 

60/201,679 and 60/222,871 filed on April 28, 2000 and August 3, 2000 respectively, 
which are both incorporated herein in their entirety. 

P. A rKCrR m TND OF INVENTION 
Computerized tomography has been used for many years, especially in 
10 the area of medical imaging, for providing nondestructive or noninvasive techniques 
for generating sectional views of an object. A good description of such systems 
appears in U.S. Patent No. 5,414,622, which is incorporated herein by reference. 
These imaging systems measure the density or metabolic activity of a section of the 
observed object (such as a patient's body) and produce raw acquisition data. This 
15 acquisition data consists of tomographic projections or is processed to obtain 
tomographic projections, which may be used in a tomographic reconstruction 
procedure to obtain image data associated with a sectional view of the observed 
object. These techniques are presently employed in medical imaging systems such as 
X-ray computerized tomography (CT) systems, positron emission tomography (PET) 
20 systems and Single Positron Emission Computerized Tomography (SPECT) systems. 
Applications in fields other than medical imaging include seismology and geophysical 
sciences where ground penetrating acoustic radiation is employed and reflected 
signals are processed to obtain tomographic projections of sub-surface structures. 
Other applications will be apparent to one skilled in the art. 
25 These systems can be described mathematically as follows. A section 

of the object to be observed by the tomographic device is represented by an image/ c 
to , x 2 ), where xj and x 2 are the spatial dimensions of the image. Jri the discrete case,/ 
is an image matrix with Ni x N 2 elements and each member of/includes pixel 
information such as, grayscale value or brightness. An estimation of/is derived by a 
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tomographic reconstruction procedure from the tomographic projections off denoted 
Y[*,a], and defined as: 

Y[t,a ] = Rf[t,a ] + W[t,a] (Eq. 1) 

where W is an additive noise, R is the discrete Radon transform, which 
5 models the tomographic projection process, a is the angle of projection and t is the 
position on the line integral. The discrete Radon transform is derived from its 
continuous version Rc, which is equivalent to the X-ray transform in two dimensions 
described in S.R. Deans, The Radon Transform And Some Of Its Applications (John 
Wiley and Sons) (1983) and expressed as 

CO 00 

10 (R c f c )( t - a )= ljf c (x„x 2 )d(x i cosa + x 2 sina-t)dx 1 dx 2 (Eq.2) 

—GO— 00 

where f c is a square integrable function 5 is the Dirac mass function 
where 5(x)=l when x=0 and S(x)=0 otherwise and a is an angle from 0 to 2n radians. 
Figure 5 is a diagram showing a typical relationship between R,/, t and a. 

One skilled in the art will recognize that there are several ways to 
1 5 define the discrete Radon transform based on the continuous Radon transform. Some 
of those methods are described in P. Toft, TJte Radon Transform and Implementation, 
Ph.D. thesis, Department of Mathematical Modeling, Technical University of 
Denmark (1996). Typically, a line integral along xj cos a + x 2 sin a = tis 
approximated by a summation of the pixel values inside the strip t-1/2 < nj cos a + n 2 
20 sin a<t+ 1/2. 

The noise W from Eq. 1 is traditionally modeled as Gaussian or 
Poisson noise. However, since the tomographic projections Y are often processed to 
incorporate various corrections, such as attenuation correction, scatter correction, 
resolution correction and geometric correction, the resulting noise W does not actually 
25 comply with such statistical models. 

Once tomographic projections are obtained, reconstruction of the 
desired image data proceeds in two steps. The steps may be performed in either order. 
The first of these steps is backprojection. This step entails computing the inverse 
radon transform (R 1 ) of the tomographic projections Y. A discrete backprojection 
30 may be directly computed by (a) radial interpolation, using techniques well known to 
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one skilled in the art and (b) deconvolution to amplify the high frequency components 
of the 1 -dimensional Fourier transforms of Y in the. direction of t by multiplying the 
1 -dimensional Fourier transforms of Y in the direction of t by a high-pass filter, such 
as a ramp filter. 

5 Second, regularization if performed. It follows from Eq. 1 that 

computing the inverse Radon transform of the tomographic projections yields: 

R- l (Y) = f + R- l (W)> (Eq.3) 
Where the desired image /is contaminated by additive noise 
Because the Radon transform R is a smoothing transform, computing its inverse R" 1 in 

10 the presence of additive noise W is an ill-posed inverse problem, so R~ 1 (W) represents 
a large additive noise. Consequently, the regularization step is necessary to remove as 
much of the additive noise as possible from the final reconstructed data. 

The prior art teaches two approaches to regularization, both of which 
suffer shortcomings. The most popular regularization technique in the prior art is 

15 known as Filtered Back-Projection (FBP). FBP and its derivatives are linear filtering 
techniques in the Fourier space. FBP suffers from limitations due to the fact that the 
vectors of the Fourier basis are not adapted to represent spatially inhomogeneous data 
typically found in images. The second well-known technique for regularization 
utilize iterative statistical models (ISM) designed to implement Expectation- 

20 . Maximum (EM) and Maximum A Posteriori (MAP) estimators. These techniques are 
described in G. McLachlan and T. Krishnan, The EM Algorithm and Extensions 
(New York: Wiley) (1997) and L.A. Shepp and Y. Vardi, Maximum Likelihood 
Reconstruction For Emission Tomography, Transactions on Medical Imaging, 1982, 
at 1 13-122. Although ISM methods can show an improvement over FBP techniques, 

25 they suffer from several drawbacks, including long computation times, lack of 

theoretical understanding to characterize estimation error and convergence problems. 

SUMMARY OF THE INVENTION 

It is an object of the present invention to provide improved methods for 
regularization of tomographic projection data. . 
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In one embodiment of the present invention, a method is provided for 
reconstructing image data corresponding to a cross-section of an object, such as a 
patient. Tomographic projection data supplied by a tomographic device is 
backprojected. The backprojected data is decomposed in a wavelet or wavelet packet 
5 vector family. The decomposed data is then diagonally operated upon and the 

modified data is then recomposed to form image data corresponding to the desired 
image by inverting the wavelet vector decomposition. In a highly preferred 
embodiment, the diagonal operation consists of thresholding the decomposed data to 
attenuate or eliminate certain decomposition coefficients. 

10 In another preferred embodiment of the present invention, the supplied 

tomographic projection data is first decomposed in a wavelet or wavelet packet vector 
family. The decomposed data is then diagonally operated upon. Next, the data is 
recomposed by inverting the wavelet vector decomposition. Finally, the data is 
backprojected to form image data corresponding to the desired image. 

15 In another embodiment, the supplied tomographic data is decomposed 

and diagonally operated upon as described above. The data is then selectively 
amplified to approximate the deconvolution step of backprojecting the data. The data 
is then recomposed by inverting the wavelet vector decomposition and the resulting 
data is radially interpolated to form image data corresponding to the desired image. 

20 In a further embodiment, the selective amplification can occur before the diagonally 
operating step. 

In a highly preferred embodiment, the wavelet or wavelet packet 
vector family chosen to decompose the tomographic data is a wavelet packet family. 
A dictionary of wavelet packet decompositions to choose is provided. One of those 

25 decompositions is selected by calculating an estimation of final estimation error based 
on the selected wavelet packet decomposition and a model of the noise in the . 
tomographic projection data for each possible wavelet packet decomposition and 
choosing the wavelet packet decomposition with the minimum estimated estimation 
error. In another preferred embodiment, the model of the noise in the tomographic 

30 projection data is extracted from the tomographic projection data. In yet another 
preferred embodiment, prior information about the desired image gained from 
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analysis of previous similar images, e.g. of phantoms", is used to calculate the 
estimated estimation error. 

In another preferred embodiment, tomographic information 
corresponding to multiple cross-sections of an object, such as a patient, is provided 
5 and the desired image is reconstructed using all of the supplied data. 

BRIEF DESCRIPTION OF THE DRAWINGS 

Figure 1 shows a flow diagram of one embodiment of a method 
according to the present invention. 

Figure 2 shows a flow diagram of another embodiment of a method 
1 0 according to the present invention. 

Figure 3 shows a flow diagram of another embodiment of a method 
according to the present invention. 

Figure 4 shows a flow diagram of a preferred wavelet family selection 
step useful in the methods of Figures 1-3. 
15 Figure 5 shows a conceptual relationship between an object to be 

imaged and a portion of its Radon transform. 

DFTATTFD DESCRIPTION OF THE PREFERRED EMBODIEMENTS 

Figure 1 shows one method of implementing the present invention. In 
step 101, tomographic data generated by a tomographic device such as an X-ray CT 

20 scanner, PET or SPECT device is supplied and backprojected. The tomographic data 
may be raw, or may be pre-processed to correct for attenuation, scatter, resolution, 
geometric anomalies, etc. These pre-processing methods are well known to one 
skilled in the art. The backprojection step, as previously described, involves radially 
interpolating the data as well as deconvolving the data. Radial interpolation is well 

25 known to one skilled in the art. In the preferred embodiment, linear interpolation is 
performed, however other radial interpolation techniques such as the nearest neighbor 
or spline radial interpolation algorithms can be utilized. Deconvolution is also well 
known to one skilled in the art and consists of multiplying the one dimensional 
Fourier transform of the projection data in the direction of t by multiplication by a 

30 high-pass filter such as a ramp filter. The entire backprojecting step is described in 
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detail in U.S. Patent Nos. 5,778,038 and 5,414,622, incorporated herein by reference. 
Attached as Appendix A is a MATLAB™ code listing implementing one embodiment 
of the present invention. The backprojecting steps are performed during the function 
filtback in the code. 

5 In step 103, a wavelet or wavelet packet vector family is selected to 

decompose the backprojected data. Either one-dimensional or two-dimensional 
wavelet vector families can be used. The wavelet or wavelet packet vector family can 
be an orthogonal basis, a biorthogonal basis or a complex or non-complex frame. A 
biorthogonal basis differs from an orthogonal basis in that the vector family used to 

1 0 reconstruct the data (denoted g* herein) is different than the vector family used to 
decompose the data (denoted g herein). Biorthogonal bases are preferable when 
particular properties of the wavelet vectors or wavelet packet vectors, such as support 
size, regularity, symmetry and number of vanishing moments are difficult or 
impossible to achieve using a orthogonal basis. Selection of a wavelet packet vector 

15 family allows a more flexible and accurate segmentation of the frequency domain, as 
the supports of the Fourier transforms of the wavelet packets nearly segment the 
frequency space. Complex wavelet frames and complex wavelet packets are 
translation invariant transforms that generate complex coefficients. These vector 
families provide improved directional selectivity. All of these vector families will be 

20 known to one skilled in the art and are described in S. Mallat, A Wavelet Tour of 
Signal Processing (Academic Press, 2d Ed.) 1999. The chosen vector family is 
denoted B herein. In a preferred embodiment demonstrated by the code in Appendix 
A, a two dimensional wavelet packet decomposition based on the Symmlet wavelet is 
used. 

25 Referring to Figure 4, a highly preferred embodiment of the vector 

family selection step 1 03, where the vector family chosen is a one or two dimensional 
wavelet packet frame, will be described. The function SPECTMacrTomoFiltback in 
the MATLAB™ code listing contained in Appendix A implements the steps involved 
in the embodiment illustrated in Figure 4 using a two dimensional wavelet packet 

30 frame and provides further description of the steps illustrated in Figure 4. In step 401, 
a model of the noise in the tomographic projection data is derived. The model may be 
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predetermined, or it may be determined from information extracted from the 
tomographic projection data. A dictionary of potential wavelet packet 
decompositions is also supplied. Generating the dictionary of wavelet packet 
decompositions is well known to one skilled in the art. In step 403, a possible wavelet 
5 packet decomposition is selected from the supplied dictionary. In step 405, an 
estimation of the final estimation error of the reconstructed image based on the 
selected wavelet packet decomposition and noise model is calculated. In step 407, the 
process repeats from step 403 if other potential wavelet packet decompositions in the 
supplied dictionary remain to be analyzed. 
10 In step 409, the wavelet packet decomposition that minimizes 

estimated final estimation error is chosen. 

An algorithm for selecting the wavelet packet which minimizes the 
estimated final error as shown in steps 405 and 409 is well known in the art and is 
described in R. Coifrnan and M. Wickerhauser, Entropy-Based Algorithms for Best 
15 Basis Selection, 38 IEEE Trans. Info. Theory, 713 (Mar. 1992) and in U.S Patent Nos. 
5,526,299 and 5,384,725 which are incorporated by reference herein. Generating a 
criterion to employ in implementing the algorithm of Coifrnan and Wickerhauser that 
estimates the estimation error based on a model on the noise in the tomographic 
projection data is well known to one skilled in the art and is described in S. Mallat, A 
20 Wavelet Tour of Signal Processing (Academic Press 2d Ed.) (1 999). Although not 
shown in Figure 4, as noted in the code listing of Appendix A the criterion can be 
based on information about the image to be observed in addition to information about 
the noise in the tomographic projection data by making reference to a prior 
observation, e.g. of a "phantom", of an object similar to the object under 
25 investigation. 

Returning to Figure 1, in step 105, the backprojected tomographic data 
is decomposed into vector family B. The algorithms employed to decompose the data 
are well known to one skilled in the art and will depend on the vector family chosen. 
In the preferred embodiment implemented in the code fisting of Appendix A, a fast 
30 filter-bank transform is used. A fast fining scheme may also be employed. Fast filter 
bank transforms are well known to one skilled in the art and are described in S. 
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Mallat, supra. Lifting schemes are also well known to one skilled in the art and are 
described in W. Sweldens, The lifting Scheme: A Construction of Second Generation 
Wavelets, 29 SIAM J. of Math. Analysis 511, 546 (1997). Techniques for 
decomposing in other wavelet vector families, using other filter bank techniques such 
as translation averaging are also possible. 

In step 107, the decomposed data is diagonally operated on. This step 
consists of multiplying the decomposed data coefficient by a scalar value pj. In a 
preferred embodiment, p s is a thresholding operator T;, defined: 
if\x\>T\ 



,or 



(Eq.4) 



10 



x-Ti if jc>7) 
x + Ti ifx<-Ti 
0 if\x\<T { 



>, or 



(Eq. 5) 



if\x\>T ( 

sign(x)a\x\-T i ) z/7^1-£)<H<7;. 



(Eq. 6) 



where TO\, or using a sigmoid function. Eq. 4 describes a soft- 
thresholding operator, Eq. 5 describes a hard thresholding operator, and Eq. 6 
describes a "piece-wise affine" thresholding operator. T (and in the case of piece 

1 5 wise affine thresholding, A.) is chosen to minimize a numerical measure of error, or to 
satisfy perceptual requirements according to the needs of the person or device that 
will interpret the reconstructed image. In a typical case, X is 1/2 and T is three times 
the estimated standard deviation of the noise in the tomographic data. 

Finally, in step 109, the data is recomposed using the inverse 

20 decomposition of the wavelet or wavelet packet vector family decomposition 

performed in step 105. The inverse decomposition process will be known to one 
skilled in the art and is based on the same techniques utilized in step 105. Thus, 
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where the forward transform utilizes filter bank techniques, the reverse transform uses 
the same techniques. In the preferred embodiment implemented in the code listing of 
Appendix A, an inverse transform using a filter bank is employed. The resultant data 
is image data corresponding to the desired cross-sectional image which may viewed 

5 or printed using an apparatus well known to one skilled in the art. 

A mathematical representation of the steps performed in steps 101 
through 1 09 follows. Let B = {gi ,g*J os i ^ with M >: Ni x N 2 - 1 be the vector family 
(previously described) in which the regularization of the data is computed where 
{g*i}i is the dual family of the family of vectors used during the reconstruction 

10 process. When {gi}i is an orthogonal basis, the families g and g* are identical. The 
diagonal estimator F (which is the matrix of regularized image data) of /is expressed: 

i 

Figure 2 illustrates a slightly modified version of the embodiment 
described with reference to Figure 1 . In the embodiment of Figure 2, the step of 

15 backprojecting the data 209 occurs last. The backprojection step 209 is identical to 
backprojection step 101 shown in Figure 1 in all other respects. The remaining steps 
shown in Figure 2 are also identical to the corresponding steps shown in Figure 1. 
Thus, step 201 corresponds to step 103. Step 203 corresponds to step 105. Step 205 
corresponds to step 107. Step 207 corresponds to step 109. Accordingly, unlike the 

20 embodiment illustrated in Figure 1, the regularization process performed in steps 201 
through 207 occurs before application of the inverse operator R" 1 . 

A mathematical representation of the steps performed in steps 201 
through 209 follows. Let B = {g h g*J o* isM with M >NiX N 2 -l be the vector family 
in which the regularization of the data is computed. The diagonal estimator F of/is 

25 expressed: 

F^R-^PAY'SA;)- CEq.8) 

Figure 3 illustrates another embodiment of the present invention, 
where the process of backprojecting the tomographic data is partially performed 
during the regularization process. In step 301, a wavelet packet vector family is 
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selected for decomposition of supplied tomographic data. The process is similar to 
that described above with reference to step 103 in Figure 1, except that this 
embodiment requires the use of a wavelet packet or complex wavelet packet frame as 
the selected vector family. 
5 La step 303, the tomographic data is decomposed using a similar 

process as described with reference to step 105 in Figure 1. 

In step 305, the data is diagonally operated on in a process similar to 
that described with reference to step 107 in Figure 1. 

In step 307, the data is selectively amplified. This step consists of 
10 amplifying selected wavelet packet coefficients to approximate amplification of the 
high frequency components of the tomographic projection data. In accordance with 
this embodiment, it is unnecessary to perform a Fourier transform and its inverse to 
backproject the tomographic projection data, rather all that is required in addition to 
the amplification step is radial interpolation described below. Selection of wavelet 
15 packet coefficients to amplify in order to approximate the traditional deconvolution 
step of backprojection will be apparent to one of skill in the art. In a preferred 
embodiment, each wavelet packet vector has a frequency support which is mostly 
concentrated on a band of the Fourier domain. In other words, the supports of the 
Fourier transforms of the various wavelet packets approximately segment the Fourier 
domain. The deconvolution by a ramp filter function is then approximated by 
amplifying the coefficients in the wavelet packet space by a value that is the mean 
value of the ramp filter function on the support of the corresponding wavelet packet 
vector. 

In step 309, the data is recomposed using the inverse decomposition of 
the wavelet packet decomposition performed in step 303. The process is similar to 
that described with respect to step 109 in Figure 1. 

Finally, in step 311, the data is radially interpolated to complete 
backprojection. This step is well known to one skilled in the art and is discussed 
above with respect to step 101 . 



-10- 



WO 01/84500 



PCT/USOl/13737 



As symbolized by arrow 306, the steps of diagonally operating 305 and 

selectively amplifying data 307 can occur in either order. In other words, step 307 

may be performed before step 305. 

The embodiments described above could also be combined to improve 

5 regularization of the tomographic data. For instance, the embodiment of Figure 2 
could be combined with the embodiment of Figure 1 , beginning with step 201, 
proceeding through Figure 2 to step 209 then moving to step 103 of Figure 1 and 
proceeding to step 109. In such an instance, the same decomposition vector family B 
could be used, or two different vector families could be selected during steps 201 and 

10 103. 

Additionally, the invention is intended to include regularization of 
three-dimensional image data in addition to traditional two-dimensional cross-section 
images. In a three dimensional embodiment of the present invention, a series of 
tomographic projections of two-dimensional cross-sectional slices of the observed 
15 object are provided by a tomographic device. The complete three-dimensional data 
set is then operated upon, which has the effect of better discriminating between 
desired information and additive noise. 

An example of a three-dimensional embodiment of the present 
invention is now described with reference to Figure 1. In step 101, data 
20 corresponding to tomographic projections of a series of cross-sectional slices of the 
observed object are provided and the projections are backprojected for each cross- 
sectional slice, generating data corresponding to a series of noisy cross-sectional 
slices of the observed object The backprojection step for each cross-sectional slice is 
identical to that previously described with respect to Figure 1 in the two-dimensional 
25 embodiment. 

Although not shown in Figure 1 , it is advantageous that a two- 
dimensional regularization on each backprojected slice be performed using one of the 
embodiments previously described before proceeding to step 103. 

In step 103, a wavelet or wavelet packet vector family is chosen to 
30 decompose the three dimensional data Either one-dimensional, two-dimensional or 
three-dimensional wavelet vector families can be used. The wavelet or wavelet 
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packet vector family can be an orthogonal basis, a biorthogonal basis or a complex or 
non-complex frame. The chosen vector family is denoted B herein. 

In step 105, the data is decomposed in the selected vector family. The 
step is similar to that previously described in the two-dimensional situation. In a 
5 highly preferred embodiment, a three-dimensional dyadic wavelet basis is chosen and 
the "a trous algorithm" is used to compute the decomposition. This algorithm is 
another method of performing a translation invariant decomposition and, like 
translation averaging, is a refinement of the filter bank decompositions previously 
discussed, although it is implemented differently than the translation averaging 
1 0 technique. The algorithm is an "undecimated" modification of the filter bank 

algorithm, without subsampling between each wavelet scale. The "a trous algorithm" 
is well known to one skilled in the art and is described in Mallet supra. An 
implementation of the "a trous algorithm" in MATLAB code is presented 
Appendix B. 

15 In step 1 07, the data is diagonally operated upon using a process 

similar to that discussed in the two-dimensional embodiment. In a preferred 
embodiment, where a three-dimensional dyadic wavelet basis is chosen to decompose 
the backprojected data, the diagonal operation consists of thresholding the wavelet 
modulus. Specifically, the dyadic transform is computed with a family of discretized 

20 translations and dilations of three wavelets y 7 , and that are the partial 

derivatives of a smoothing function ((>. Smoothing functions suitable for this purpose 
are well known to one skilled in the art and are discussed in Mallet supra. A typical 
smoothing function is a quadratic spline. Thus: 

y -fc,^). flfajaai , ¥ 3 ( ,„ X! ^)=3fe 

dx } dx 2 ox 3 

25 . 

The three-dimensional vector family {g\}i is then written {\\f k j}k=i,2,3j- 
For a given j, the three wavelets have been equally dilated and are translated on the 
same position, but they respectively have a horizontal, vertical and axial direction. 

The decomposed image /has three components which can be written 
30 as frame inner products, as: 
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T*jf = {fW k j). k = l,2,3. (Eq.9) 
Because \|/, \|/ 2 and \\r 3 are partial derivatives of <|), these components 
are proportional to the coordinates of the gradient vector off smoothed by a dilated 
version of <|>. From these coordinates, a computation of the angle of the gradient 
5 vector is made, which indicates the direction in which the partial derivative of the 
smoothed image/has the largest amplitude. The amplitude of this maximum partial 
derivative is equal to the modulus of the gradient vector and is therefore proportional 
to the wavelet modulus, expressed as: 

Mjf = i\T}f\ 2 + \Tlf\ 2 +\T]f\ 2 . (Eq. 10) 

1 0 Finally, thresholding of the modulus Mj/is performed, rather than 

thresholding each wavelet transform component T k j/ This is equivalent to selecting a 
direction in which the partial derivative is a maximum and thresholding the amplitude 
of the partial derivative in this direction. This constitutes an adaptive choice of the 
wavelet direction in order to best correlate the signal. The coefficients of the dyadic 
1 5 wavelet decomposition are then computed back from the thresholded modulus and the 
angle of the gradient vector. 

In step 109 the data is recomposed to generate a series of regularized 
cross-sectional slices. 

A mathematical representation of the steps performed in steps 101 
20 through 109 for the situation with N3 number of cross-sectional slices follows. Let B 
= {gi,g*i}i be the vector family in which the regularization of the data is computed. 
The three-dimensional slice-by-slice backprojected data is: 

V0<z<Ar 3 -l X[ 9 .,z] = R- l (Y[,.,z]). (Eq. 11) 

where z corresponds to the axial direction of the object under investigation. The 
25 diagonal estimator F of /is : 

F = ^ Pi {(X, gi ))g;. (Bq.12) 

I 

It will be apparent to one skilled in the art that the embodiment 
illustrated in Figure 2 could be adapted in a similar manner to operate on three 
dimensional data as well. In that case, a mathematical .representation of the steps 
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performed in steps 201 through 209 for the situation with N 3 number of cross- 
sectional slices follows. The regularized tomographic projection data is given by: 

r = £P,-((^)k : (Eq-13) 
and the estimator F of/is: 
5 \/0<z<N 3 -l F = [.,.. z]=R- l ^Y[,.,z]j. (Eq. 14) 

The foregoing merely illustrates the exemplary embodiments of the 
invention. Various modifications and alterations to the described embodiments will 
be apparent to those skilled in the art in view of the teachings herein. It will thus be 
fully appreciated that those skilled in the art will be able to devise numerous systems 
10 and methods which, although not explicitly shown or described, embody the 

principles of the invention and are thus within the spirit and scope of the invention, 
the scope of the invention being defined only by the appended claims. 



01B4S00A2J_> 
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APPENDIX A 

function x = Af fine Thresh (y, t) 
% AffineThresh Apply Affine Threshold 
% Us age 
5 % x = Aff ineThresh(y, t) 
% Inputs 

% y Coeff before thresholding 

% t Threshold 
% Outputs 
10 % x Thresholded coeff 

% 



15 



x = sign (y) . * (abs (y) > t/2) . * (abs (y) <= t) . * (2*abs (y) -t ) + 
y .* (abs (y) > t) ; 

% J.Kalifa, 2000 



function [basis , value] = Best2dBasis (sqtree , D) 

% Best2dBasis -- Coif man- Wicker hauser Best-2d-Basis Algorithm 

20 % Usage 

% [btree, vtree] = Best2dBasis (sqtree, D) 

% Inputs 

% sqtree stat-quadtree (output by Calc2dStatTree) 

% D maximum depth of tree-search 

25 % Outputs 

% btree basis-quadtree of best basis 

% vtree value of components of best basis 

% vtree (1) holds value of best basis 

% 

30 % Description 

% The best-basis algorithm is used to pick out the ""best 1 1 

% basis from all the possible bases in the packet table. 

% Here ""best 1 ' means minimizing an additive measure of 

% information, called entropy by Coif man and Wickerhauser . 

35 % 

% Once the stat-quadtree of entropy values is created, Best2dBasis 

% selects the best basis using the pruning algorithm described in 

% Wickerhauser 1 s book. 

40 % Examples 

% % get best basis for "image 1 ' 

% % use Coiflet 3 filter and CW Entropy 

% [n,D] = quadlength (image) ; 

% qmf = MakeONFilter ( "Coiflet 1 , 3) ; 

45 % sqtree = Calc2dStatTree ( 1 WP 1 , image , D , qmf , 1 Entropy 1 ) ; 

% [btree, vtree] = Best2dBasis (sqtree, D) ; 

% 

% See Also 

% Calc2dStatTree, FPT2 CP, FPT2 WP 

50 % " 

% References 

% Wickerhauser, M.V. __Adapted_Wavelet_Analysis_. AK Peters 

(1994) . 

% 



55 



global W1.VERBOSE 
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scan 



basis = zeros (size(sqtree) ) ; 
value = sqtree; 
for d=D-l:-l:0, 

for bx=0 : (2 A d-l) , 

for by=0 : (2^d-l) , % scan tree, bottom->top left->right 



10 



15 children, RU486 



20 



vparent = sqtree (qnode (d,bx,by) ) ; 

vchild = value (qnode (d+1, 2 *bx ,2*by)) + ... 

value (qnode (d+1, 2*bx+l,2*by) ) + ... 

value (qnode (d+1, 2 *bx , 2*by+l)) + 
value (qnode (d+1, 2*bx+l, 2*by+l) ) ; 

if (vparent <= vchild) , % if parent better than 



else 



end 



basis (qnode (d,bx, by) ) = 0; 

value (qnode (d, bx, by) ) = vparent ; 

basis (qnode (d,bx, by) ) = 1; 

value (qnode (d, bx, by) ) = vchild; 



25 



end 

end 
end 

if strctnp (WL VERBOSE, 'Yes ' ) 

fprintf ( 'best basis %g \n », value (1) ) 
end 



30 % 

% Copyright (c) 1993. Jonathan Buckheit and David Donoho 
% 



35 % 

% Part of WaveLab Version 8 02 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 
40 % 



function s = cyclespin2 (x, i , j ) 

[11, 12] = size (x) ; 
45 z = x( (11+1-i) :11, :) ; 

z(i+l:ll,:) = X(l: (11-i) , : ) / 

S = Z ( : , (12 + 1-j ) :12) ; 

S ( : , j+l:l2) = Z(: ,1: (12-j) ) ; 

50 % Written by Maureen Clerc and Jerome Kalifa, 1997 

% clerc®cmapx.polytechnique . f r , kalif a@cmapx .polytechnique . f r 



55 function [maxi , tree , T] = Dis2dStatTree (img,D,qmf , I , mult_sigma) 
[n, J] = quadlength(img) ; 
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boxlen = n; 
boxcnt = 1; 
coef = img; 
coefl - I; 

for d = 0 :D, 

for bx=0 : (boxcnt- 1) , 
for by=0: (boxcnt- 1) , 



10 [lox, nix, loy, hiy] « quadbounds (d, bx,by, n) ; 

Quad = coef (lox: hix, loy : hiy) ; 
QuadI = coefl (lox: hix, loy: hiy) ; 



15 



[si, s2] =size (Quad) ; 



T_node=sqrt (4*log (n) ) *sqrt (1/ (sl*s2)*sum (sum (QuadI . *QuadI) ) ) ; 
T_node=mult_sigma*sqrt (1/ (sl*s2) * sum (sum (QuadI . *QuadI) ) ) ; 
rl= abs(Quad) < T_node; 
r2=l-rl; 
20 pl=Quad."2; 

ptl=pl.*rl; 
p2=QuadI. A 2; 
pt2=p2 . *r2; 

tree (qnode (d, bx, by) ) = sum(sum(ptl) ) +sum (sum (pt2) ) ; 
25 T ( qnode ( d , bx , by) ) = Tjatode; 

maxi (qnode (d, bx, by) ) =max (max (Quad) ) ; 



% fin du bloc 

30 

if d < D, % prepare for finer scales 
quad = DownQuad (Quad , qmf , by , bx) ; 
quadl = DownQuad (QuadI , qmf , by , bx) ; 
coef (lox: hix, loy : hiy) = quad; 
35 coefl (lox : hix, loy : hiy) = quadl; 

end 
end 
end 

boxlen = boxlen/2; 
40 boxcnt = boxcnt *2 ; 

end 



function out = DiscrDenoise2 Radon (basis , img, T, qmf , type) 
% DiscrDenoise2 -- Denoising of inverse filtered image 
45 % in best discriminant basis 

% Usage 

% out = DiscrDenoise2 (basis, img, T, qmf ,D, maxi) 
% 

% See Also 
50 % WPAnalysis, DisBasis, Discr 
% 

[n,J] = quadlength(img) ,- 
coef = img; 

% 

55 % initialize tree traversal stack 
% 

stack = zeros (3,4* (J+l) +1) ; 
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k = 1; 

stack (:,k) = [0 0 0]'; % d, bx, by 

% 

while (k > 0) , 

5 

% pop stack 
d = stack (1, k) ; 
bx = stack (2 ,k) ; 
by = stack ( 3, k) ; 
10 k-k-1; 

[lox hix loy hiy] = quadbounds (d, bx, by , n) ; 
Quad = coef (lox: hix, loy :hiy) ; 

if (basis (qnode (d,bx, by) ) — 0) , % nonterminal node 
15 quad = DownQuad (Quad,qmf ,bx,by) ; 

coef (lox: hix, loy : hiy) = quad; 

% push stack 

k = k+1; stack(:,k) = [(d+1) (2*bx) (2*by) ]■; 
20 k = k+l; stack(:,k) = [(d+l) (2*bx+l) (2*by) ]»; 

k « k+1; stack(:,k) = [(d+l) <2*bx) (2*by+l) ]•; 
k « k+1; stacker ,k) = [(d+l) (2*bx+l) (2*by+l) ]»; 

else 

coef (lox:hix, loy :hiy) =Quad; 
25 thresh=T (qnode (d,bx,by) ) ; 

if strcmp (type, 1 A 1 ) , 

coef (lox: hix, loy : hiy) = Af f ineThresh (Quad, 

thresh) ; 

end 

30 if strcmp (type, 'S' ) , 

coef (lox:hix, loy:hiy) « Sof tThresh (Quad, thresh); 
end 

if strcmp (type, »H' ) , 

coef (lox: hix, loy : hiy) = HardThresh (Quad, thresh); 
35 end 
end 

end 

out = IPT2_WP (basis, coef ,qmf) ; 



40 



function d = DownDyadHi (x, qmf ) 

% DownDyadHi Hi -Pass Downsampling operator (periodized) 
% Usage 

% d = DownDyadHi (x, f) 
45 % Inputs 

% x 1-d signal at fine scale 
% f filter 
% Outputs 

% y 1-d signal at coarse scale 



50 



% See Also 

% DownDyadLo, UpDyadHi , UpDyaciLo, FWT_PO, iconv 
% 



d = iconv ( MirrorFilt (qmf) , lshif t (x) ) ; 
55 n = length (d) ; 

d - d(l:2: (n-l) ) ; 
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% 

% Copyright (c) 1993. Iain M. Johnstone 
% 



% 

% Part of WaveLab Version 802 

% This is Copyrighted Material 

% For Copying permissions see COPYING. m 

10 % Comments? e-mail wavelab@stat.stanford.edu 
% 



function d = DownDyadLo (x, qmf ) 
15 % DownDyadLo -- Lo-Pass Downsampling operator (periodized) 
% Usage 

% d = DownDyadLo (x, f) 
% Inputs 

% x 1-d signal at fine scale 
20 % f filter 
% Outputs 

% y 1-d signal at coarse scale 
% 

% See Also 

25 % DownDyadHi, UpDyadHi , UpDyadLo, FWT_P0, aconv 



30 



% 



d = aconv (qmf ,x) ; 

n = length (d) ; 

d = d(l:2 : (n-1) ) ; 



% Copyright (c) 1993. Iain M. Johnstone 



35 

% 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
40 % Comments? e-mail wavelab@stat.stanford.edu 
% 



function q = DownQuad (Quad, qmf ,xbit ,ybit) 

45 % Down Quad Split 2 -d image into 4 subbands 

% Usage 

% q = DownQuad (Quad, qmf ,xbit,ybit) 

% Inputs 

% Quad Quadlet to be split into four subbands 

50 % qmf orthogonal quadrature mirror filter 

% xbit x-position of quadlet in 2-d freq plane 

% ybit y-position of quadlet in 2-d freq plane 

% Outputs 

% q quadlet after splitting 
55 % 

% Description 

% A 2-d signal is split into four subbands: (Lo_x,Lo_y), 
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% ( Lo_x , Hi_y ) , (Hi_x , Lo_y ) , (Hi__x , Hi_y ) . 
% 

% This routine is called by other, higher-level WaveLab functions 
% and may not be useful for many users . 
5 % 

% See Also 

% FPT2_WP, UpQuad 

% 

10 [nr,nc] = size (Quad) ; 

q = zeros (nr,nc) ; 
ybit = rem(ybit,2); 
xbit = rem (xbit, 2); 

% 

15 bot = 1: (nc/2) ; 

top = bot + (nc/2) ; 

% 

for i=l :nr, 

lpchan = DownDyadLo (Quad (i , :) , qmf ) ; 
20 hp chan = DownDyadHi (Quad (i , : ) , qmf) ; 

% lpchan = ddc (Quad(i, : ) , qmf ) ; 

% hpchan = ddd (Quad (i , : ) , qmf) ; 

if xbit, 

q(i,top) = lpchan; 
25 q(i,bot) = hpchan; 

else 

q(i,top) = hpchan; 
q(i,bot) = lpchan; 

end 

30 end 
% 

bot = 1: (nr/2) ; 
top = bot + (nr/2) ; 

% 

35 for j=l:nc, 

lpchan = DownDyadLo (q ( : , j ) 1 , qmf) 1 ; 
hpchan = DownDyadHi (q ( : , j ) ' , qmf) 1 ; 
% lpchan = ddc (q ( : , j ) 1 , qmf) 1 ; 

% hpchan = ddd (q ( : , j ) ' , qmf) 1 ; 

40 if ybit, 

q(top,j) = lpchan; 
q(bot,j) = hpchan; 

else 

q(top,j) = hpchan; 
45 q(bot,j) = lpchan; 

end 

end 

50 % 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 
55 % 
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function 

recim=filtback (sinogram, filter type ,X, deltarho , del tax, rhowindow, interp 
ol) 

% filtback 
5 % 

% A routine for computing the inverse Radon transform from a 
(rho, theta) 
% sinogram. 
% 

10 % Peter Toft 

% Department of Mathematical Modelling, 
% Technical University of Denmark 
% 1997 

15 global INTERPOL FILTERTYPE DEBUGLEVEL 

if length (DEBUGLEVEL) ==0 

DEBUGLEVEL=0; 
end 

20 

if nargin==0 

sinogram= zeros (101, 100) / 

sinogram ( (1+size (sinogram,!) ) /2, :) =ones (1, size (sinogram, 2) ) ; 
here=0 
25 if DEBUGLEVEL>0, 

figure; 

image sc (sinogram) ; 
colorbar ; 

colormap (gray (64) ) ; 
30 title ( f Demo sinogram 1 ); 

drawnow 
end 
end 

35 [NoRho, NoTheta] =size (sinogram) ; 

if nargin<2 

if length (FILTERTYPE) >0 , 
f iltertype=FILTERTYPE ; 
40 if DEBUGLEVEL >0 

disp( 1 Using global setting of the variable f iltertype ' ) ; 

end 
else 

f iltertype= ' ramp ' ; 
45 if DEBUGLEVEL >0 

fprintf ( 1 Using default setting of the variable f iltertype 
%s\n» , f iltertype) ; 
end 
end 
50 end 

if nargin<3 
X=NoRho ; 
if DEBUGLEVEL >0 

55 fprintf ( 1 Using default setting of the variable X = %i\n , ,X); 

end 
end 
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if nargin<4 
deltarho=l; 
if DEBUGLEVEL > 0 

5 fprintf ( 'Using default setting of the variable deltarho = 

%f\n' , deltarho) ; 
end 

end 

10 if nargin<5 

deltax= deltarho ; 
if DEBUGLEVEL >0 

fprintf ( 'Using default setting of the variable deltax = 
%f \n' , deltax) ; 
1 5 end 
end 

if nargin<6 
rhowindow=0 / 
20 if DEBUGLEVEL >0 

fprintf ( 'Using default setting of the variable rhowindow = 
%i\n' , rhowindow) ; 

end 
end 

25 

if nargin<7 

if 1 eng t h ( INTERPOL ) > 0 , 
interpol=INTERPOL ; 
if DEBUGLEVEL > 0 

30 disp ('Using global setting of the variable interpol ') ; 

end 
else 

interpol =1 ; 

if DEBUGLEVEL >0 

35 fprintf ( 'Using default setting of the variable interpol = 

%i\n ' , interpol) ; 
end 
end 
end 



40 



% Harm windowing to avoid edge effects 



i f rho windows = 1 , 

HannWindow=0 .5* (1-cos (2*pi/ (NoRho-1) * (0 :NoRho-l) ) ) ' ; 
45 sinogram= sinogram. * (HannWindow*ones (l r NoTheta) ) / 

end 

% Filtering 

50 %Valid for nearest neighbour and linear interpolation only 
maxrho=sqrt (2.0)* (X-l) /2*deltax; 

NoRhoExt=max (1+2* ceil (maxrho/deltarho+1) ,NoRho) ; 
55 half addlNoRho=round (0 . 5* (NoRhoExt -NoRho) ) ; 
sinogramext = zeros (NoRhoExt , NoTheta) / 

-22- 



RM.Qnmirv <wn madsonA? i > 



PCT/US01/13737 

WO 01/84500 



sinogramext (1+half addlNoRho :NoRho+half addlNoRho, : ) ^sinogram; 

NoRhoExth= (NoRhoExt-1) /2; 
if strcmp ( lower (filtertype) , 'ramp , )==l/ 
5 disp( 'Using Ramp filter 1 ) 

f ilterf- [0 -.NoRhoExth] *0 . 5/NoRhoExth/deltarho; 
elseif strcmp (lower (filtertype) , ' shepplogan ' ) ==1, 
disp( 'Using Shepp Logan filter'); ^ 
f ilterf =sin (pi/2* [0 :NoRhoExth] /NoRhoExth) /pi/deltarno; 
10 elseif strcmp (lower (filtertype) , 'sine')==l, 
disp( 'Using sin filter')/ 

f ilterf -sin (pi* [0 :NoRhoExth] /NoRhoExth) /pi/deltarho*0 . 5 ; 
else 

disp ( 'Unknown method. ' ) 
15 end 

f ilterf unc= [f ilterf f ilterf (NoRhoExth : -1 : 1) ] ' ; 

20 fSnogramLt=real(ifft(fft(s . M f ilterf unc*ones (1 ^oTheta) ) 

)); 

%if DEBUGLEVEL> 1 , 
% figure 
25 % imagesc (f sinogramext) 
% colorraap gray 
% colorbar 
% drawnow ; 
%end 



30 



45 



% Backpro j ection 
recim=zeros (X,X) ; 



relxar= ( (0 :X-1) - (X-l) /2) *deltax/deltarho; 
35 thetaindices= (0 :NoTheta-l) ; 

thetaar=thetaindices/NoTheta*2*pi; 

cosar=cos (thetaar) ; 
sinar=sin(thetaar) ; 
xcosar=cosar 1 *relxar; 
40 ysinar=sinar ' *relxar+ (NoRhoExt-1) /2 ; 
thetaimat=thetaindices ' *ones (1,X) ; 



%xcosar (t,m) =x(m) *cos (theta(t) ) 
%ycosar(t,m) =y(m) *sin (theta (t) ) +of f set 



if interpol==0, 
for xi=l:X, 

% fprintf ('%4i (%i) \r ' , xi ,X) ; 
rhomat=ysinar+xcosar ( : ,xi) *ones (1,X) ; 
50 rhoimat=round(rhomat) ; 

indmat = rho ima t + 1 +NoRhoExt * t he t a ima t ; 
for yi=l:X, , 

recim(yi,xi) =sum(f sinogramext (mdmat ( : ,yi) ) ) ; 

end 
55 end 
else 

for xi=l:X, 
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% fprintf( ! %4i (%i) \r 1 ,xi,X) ; 

rhomat=ysinar+xcosar { : ,xi) *ones (1,X) ; 
rhoimat=f loor (rhomat) ; 
drho:=rhomat-x~hoimat ; 
5 indmat=rhoimat+l+NoRhoExt *thetaimat ; 

for yi=l:X, 

recim (yi,xi) =sum(drho ( : , yi) . *f sinogramext (1+indmat ( : ,yi) ) + (1 
drho (:,yi)) . *£ sinogramext (indmat ( : ,yi) ) ) ; 
end 
10 end 
end 

recim=recim*pi/NoTheta; 

15 recim ( : , : ) =recim (X : -1 : 1 , : ) ; 

f print f ( 'Done with filtered backproj ection\n » ) ; 

if DEBUGI J EVEI J >1 / 
20 figure; 

imagesc (recim) ; 
colorbar; 

colormap (gray (64) ) / 
title ( * Reconstructed image 1 ) ; 
25 drawnow; 
end 



function x = HardThresh (y, t) 
% HardThresh Apply Hard Threshold 
30 % Usage 

% x = HardThresh (y,t) 
% Inputs 

% y Noisy Data 

% t Threshold 
35 % Outputs 

% x y l_{|y|>t} 

% 

x = y . * (abs (y) > t) / 

40 % 

% Copyright (c) 1993-5. Jonathan Buckheit, David Donoho and Iain 

Johnstone 

% 

45 

% 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
50 % Comments? e-mail wavelab@stat.stanford.edu 
% 



function y = iconv(f,x) 
55 % iconv Convolution Tool for Two-Scale Transform 
% Usage 

% y = iconv (f,x) 
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% Inputs 
% f filter 
% x 1-d signal 
% Outputs 
5 % y filtered result 
% 

% Description 

% Filtering by periodic convolution of x with f 
% 

10 % See Also 

% aconv, UpDyadHi, UpDyadLo, DownDyadHi , DownDyadljo 
% 

n = length (x) ; 
p = length (f) ; 
15 if p <= n, 

xpadded = [x ( (n+l-p) :n) x] ; 
else 

z = zeros (1 ,p) ; 

for i=l:p, 

20 imod - 1 + rem(p*n -p + i-l,n); 

z (i) = x(imod) ; 

end 

xpadded = [z x] ; 

end 

25 ypadded = filter (f , 1 , xpadded) ; 

y = ypadded ( (p+1) : (n+p) ) ; 

% 

% Copyright (c) 1993. David L. Donoho 
30 % 



% 

% Part of WaveLab Version 802 

35 % This is Copyrighted Material 

% For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 

% 



40 

function img = IPT2_WP (basis , coef , qmf) 

% IPT2_VfP -- Synthesize image from 2-d wavelet packet coefficients 



% 


Usage 




% 


img = IPT2_WP (btree , coef , qmf ) 


% 


Inputs 




% 


btree 


Quad Tree indicating wavelet packet to use 


% 


coef 


2-d wavelet packet coeffts in given basis 


% 


qmf 


quadrature mirror filter 


% 


Outputs 




% 


img 


2-d image whose 2-d wavelet packet coefft r i 


% 




basis are coef 


% 


Description 


% 
% 


Perform 


the inverse operation of FPT2_WP. 



55 [n,J] = quadlength(coef ) ; 

img = coef; 
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% initialize tree traversal stack 
stack = zeros (4,4* (J+D +1) ; 
k = 1; 

5 stack (:,k) = [0 0 0 0] » ; % d, bx, by, unmarked 

% 

while (k > 0) , 

% pop stack 
10 d = stack (l,k) ; 

bx = stack(2,k); by = stack(3,k); 

marked = stack(4,k); k=k-l; 

% fprintf(«d bx by 1 ); disp([d bx by]) 

15 if (basis (qnode (d,bx,by) ) -= 0) , % nonterminal node 

if (marked == 0) , 

% first visit, because unmarked 
20 % pushdown (marked) self 

% pushdown unmarked children 



25 0] 
o] 



30 



40 



45 



(2*by+l) 0] 
(2*by+l) 0] 



k = 
k = 


k+1; 
k+1; 


stack ( : 
stack ( : 


,k) 
,k) 


= [d bx by 1] 1 ; 
= [(d+1) (2*bx) 


(2*by) 


k = 


k+1; 


stack ( : 


,k) 


= [(d+1) 


(2*bx+l) 


(2*by) 


k = 


k+1; 


stack ( : 


,W 


- [(d+1) 


(2*bx) 




k = 


k+1; 


stack ( : 




= [(d+1) 


(2*bx+l) 





else 



% second (&last) visit, because marked 
35 % reconstruct apply UpQuad operator 



loy hiy] ) 



[lox hix loy hiy] = quadbounds (d, bx, by, n) ; 
%fprintf ( • [lox hix loy hiy] » ) ; disp ( [lox hix 



quad = img(lox:hix, loy :hiy) ; 
Quad - OpQuad(quad,qmf ,by,bx) ; 
img (lox : hix, loy :hiy) = Quad; 



end 

end 



end 



% 

50 % Copyright (c) 1993. David L. Donoho 
% 



% 

55 % Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING . m 
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function f - MakeONFi Iter (Type, Par) 
5 % MakeONFilter Generate Orthonormal QMF Filter for Wavelet 
Transform 
% Usage 

% qmf = MakeONFilter (Type, Par) 
% Inputs 

10 % Type string, 'Haar', 'Beylkin', 'Coiflet', 1 Daubechies ' , 
% 1 Symmlet 1 , ' Vaidyanathan 1 , 1 Battle 1 

% Par integer, it is a parameter related to the support and 
vanishing 

% moments of the wavelets, explained below for each 

15 wavelet. 
% 

% Outputs 

% qmf quadrature mirror filter 
% 

20 % Description 

% The Haar filter (which could be considered a Daubechies -2) was 
the 

% first wavelet, though not called as such, and is discontinuous. 
% 

25 % The Beylkin filter places roots for the frequency response 
function 

% close to the Nyquist frequency on the real axis . 
% 

% The Coiflet filters are designed to give both the mother and 
30 father 

% wavelets 2*Par vanishing moments; here Par may be one of 1,2,3,4 

or 5 . 

% 

% The Daubechies filters are minimal phase filters that generate 
35 wavelets 

% which have a minimal support for a given number of vanishing 
moments . 

% They are indexed by their length, Par, which may be one of 
% 4,6,8,10,12,14,16,18 or 20. The number of vanishing moments is 
40 par/2. 
% 

% Symmlets are also wavelets within a minimum size support for a 
given 

% number of vanishing moments, but they are as symmetrical as 
45 possible, 

% as opposed to the Daubechies filters which are highly 
asymmetrical . 

% They are indexed by Par, which specifies the number of vanishing 
% moments and is equal to half the size of the support. It ranges 
50 % from 4 to 10. 
% 

% The Vaidyanathan filter gives an exact reconstruction, but does 
not 

% satisfy any moment condition. The filter has been optimized for 
55 % speech coding. 



% 
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% The Battle -Lemarie filter generate spline orthogonal wavelet 
basis . 

% The parameter Par gives the degree of the spline. The number of 

% vanishing moments is Par+1 . 

% 

% See Also 

% FWT_PO, IWT_P0, FWT2_PO, IWT2_PO, WPAnalysis 
% 

% References 

% The books by Daubechies and Wickerhauser . 
% 

if strcmp (Type, l Haar«) , 

f = [11] . / sqrt (2) ; 

end 

if strcmp (Type, 'Beylkin 1 ) , 
f = [ .099305765374 



.424215360813 



] ; 



.449718251149 
. 026900308804 
-.088543630623 
- . 017460408696 
. 001484234782 



- .110927598348 
.155538731877 
. 019679866044 
- . 014365807969 
- . 002736031626 



699825214057 

- .264497231446 



- . 017520746267 
. 042916387274 
. 010040411845 
. 000640485329 



end 



if strcmp (Type, 1 Coif let' ) , 
if Par==l , 

f = [ .038580777748 



.126969125396 



. 077161555496 



.607491641386 
3 ; 



.745687558934 



. 226584265197 
end 

if Par==2 , 

f = [ .016387336463 



.417005184424 
. 023680171947 
] ; 



-.041464936782 ' -.0673 72 554722 
.3 8611006682 3 . 812 723 63 5450 

- . 07 64885 9907 8 - . 05943441864 6 

.005611434819 -.0 0182 32 08 871 



. 000720549445 
end 

if Par==3 , 

f = [ -.003793512864 



.405176902410 
. 071799821619 

. 015880544864 



. 007 7 82 59642 6 . 0234 52 696142 

- . 065771911281 - . 061123390003 

.793777222626 .428483476378 
-.082301927106 .034555027573 
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. 001117518771 



] ; 



- . 009007976137 
000466216960 



- . 002574517688 
- . 000070983303 



. 000034599773 
end 

if Par==4 , 

f = [ .000892313668 



.001629492013 



.007346166328 



. 081266699680 

. 782238930920 

. 096220442034 

. 015211731527 

.001266561929 
.000062339034 

] ; 



. 016068943964 
- . 056077313316 
.434386056491 
. 039334427123 
-.005658286686 
- . 000589020757 
. 000031229876 



. 026682300156 

.415308407030 

- . 066627474263 

.025082261845 

.003751436157 

- . 000259974552 

-.000003259680 



. 000001784985 
end 

if Par==5, 

f = [ -.000212080863 



.000358589677 



.002178236305 



. 023408156762 

. 052043163216 

.437991626228 
. 041289208741 

. 009164231153 

. 001662863769 

. 000140541149 
. 000003734597 



- .004159358782 
028168029062 
421566206729 
.062035963906 
032683574283 

.006764185419 

- .000638131296 

- .000041340484 
. 000002063806 



- .010131117538 
- . 091920010549 
. 774289603740 
-.105574208706 
- . 019761779012 
.002433373209 
.000302259520 
-.000021315014 
-.000000167408 



. 000000095158 
end 

end 



] ; 



if strcmp (Type, 1 Daubechies ' ) , 
if Par==4, 

f = [ .482962913145 .83 65163 0373 8 

. 224143 868042 - . 1294 09522551 

end 

if Par==6, 

f = [ .332670552950 .806891509311 

.459877502118 -.135011020010 
- .085441273882 .03522 62 91882 

end 



] ; 



] ; 
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if Par==8, 

f = [ .230377813309 
. 630880767930 
- . 187034811719 
. 032883011667 

end 

if Par==10, 

f = [ .160102397974 



.714846570553 
- . 027983769417 
.030841381836 
- . 010597401785 



] ; 



. 032244869585 
. 012580751999 



. 603 8292 69797 . 7243 0852843 8 

.138428145901 -.242294887066 

.077571493840 -.006241490213 

. 003335725285 
] ; 



end 

if Par==12, 

f = [ .111540743350 



. 129766867567 
.031582039317 



.4 94623 8 903 98 .751133 908021 

.315250351709 -.2262 64693 965 

.097501605587 .027522 86553 0 

.000553842201 .004777257511 



] ; 



. 001077301085 
end 

if Par==14, 

f = [ .077852054085 



.396539319482 



.729132090846 



.224036184994 
.038029936935 

. 000429577973 
] ; 



.469782287405 
. 071309219267 
- . 016574541631 
- . 001801640704 



- . 143906003929 
. 080612609151 
. 012550998556 
. 000353713800 



end 

if Par==16, 

f » [ .054415842243 



.312871590914 



.675630736297 



.284015542962 

.017369301002 

. 008746094047 
. 000675449406 



.585354683654 
. 000472484574 
- . 044088253931 
. 004870352993 
, 000117476784 

] ; 



- . 015829105256 
. 128747426620 
. 013981027917 
- . 000391740373 



end 

if Par==18, 

f = [ .038077947364 



.293273783279 



.243 834674613 . 604823123690 

.657288078051 .1331973 85825 
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. 030725681479 
. 022361662124 
. 001847646883 



096840783223 
067632829061 
004723204758 
000230385764 



] ; 



.148540749338 
. 000250947115 
-.004281503682 
-.000251963189 



. 000039347320 
end 

if Par==2 0, 

f = [ .026670057901 



188176800078 



.527201188932 



.249846424327 

. 093057364604 
. 033212674059 
. 001395351747 

.000116466855 

] J 

end 

end 



.688459039454 
- .195946274377 
- . 071394147166 
. 003606553567 
. 001992405295 
.000093588670 



.281172343661 
.127369340336 
-.029457536822 
-.010733175483 
- . 000685856695 
- . 000013264203 



if strcmp (Type , 1 Symmlet ' ) , 
if Par==4 , 

f = [ -.107148901418 



140317624179 



- . 041910 965125 .70373 90 6 8656 

1.136658243408 .4212345342 04 

- . 017824701442 . 04557 0345896 



3 ; 



end 

if Par==5, 

f = [ .038654795955 



. 896581648380 
. 029842499869 



. 041746864422 - . 055344186117 

.281990696 8 54 1 . 023 052 966894 

.023478923136 -.247951362613 



. 027632152958 
] ; 

end 

if Par==6 , 

f = [ .021784700327 



.004936612372 



- .166863215412 



1. 113892783926 
. 029783751299 



- . 068323121587 
.477904371333 
. 063250562660 



.694457972958 
-.102724969862 
. 002499922093 



,011031867509 
end 

if Par==7 , 



] ; 
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. 070078291222 

1. 085782709814 
.152463871896 

]; 



.043155452582 
.024665659489 
.408183939725 
. 005671342686 



.096014767936 
.758162601964 
- .198056706807 
. 014521394762 



end 

if Par==8, 

f = [ .002672793393 



- . 000428394300 



- . 021145686528 



. 038493521263 

1.099106630537 
.202648655286 
. 000766690896 



.005386388754 

- . 073462508761 

.680745347190 

. 010758611751 

- . 004783458512 
] ; 



. 069490465911 
.515398670374 
- . 086653615406 
. 044823623042 



end 

if Par==9, 

f « [ .001512487309 



.000669141509 



- . 014515578553 



. 025786445930 

. 873048407349 
.077172161097 
.016303351226 



.012528896242 
- .270893783503 
1. 015259790832 
. 000825140929 
- . 018769396836 

] ; 



. 087791251554 
. 049882830959 
.337658923602 
. 042744433602 
. 000876502539 



.001981193736 
end 

if Par==10, 

f = [ .001089170447 



,000135245020 



. 012220642630 



. 016418869426 
.667071338154 
.050256540092 

.008152816799 
. 006495728375 
1 ; 

end 

end 



- . 002072363923 
.225558972234 
088251530500 
. 045240772218 
028786231926 
. 000080661204 



. 064950924579 
- . 100240215031 
.542813011213 
.070703567550 
- . 001137535314 
- . 000649589896 
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if strcmp (Type , ' Vaidy ana than ! ) , 

f = [ -.000062906118 .000343631905 



- . 000453956620 



- . 000944897136 
- . 008839103409 
- .014853448005 
. 055892523691 
. 131971661417 
- .263494802488 
. 572797793211 



.002843834547 

. 003153847056 

- . 035470398607 

-.077709750902 

.135084227129 

.201612161775 

.250184129505 



.000708137504 

.019687215010 

.038742619293 

- .083928884366 

-.194450471766 

.635601059872 

.045799334111 



] / 



end 



if strcmp (Type , ' Battle 1 ) , 
if Par == 1, 

g = [0.578163 



0 . 012003 

0 .000652922 



0.280931 -0.0488618 -0.0367309 ... 
0.00706442 -0.00274588 -0.00155701 ... 
0.000361781 -0.000158601 -0.0000867523 



] ; 



end 



if Par == 3, 



[0.541736 0 
0.0226846 

0.00614143 0.00579932 
0.00154624 0.00133086 
0.000395946 0.000326749 
0.000103307 

] ; 



30683 -0.035498 -0.0778079 ... 

0.0297468 -0.0121455 -0.0127154 . 

-0 . 003 07 863 -0 . 0 02745 29 
-0.000780468 -0.00065562 .. 
-0.000201818 -0.000164264 . 



end 

if Par == 5, 

g = [0.528374 0.312869 -0.0261771 -0.0914068 ... 

0.0208414 0.0433544 -0.0148537 -0.0229951 

0.00990635 0.0128754 -0.00639886 -0.00746848 . 

0.00407882 0.00444002 -0.00258816 -0.00268646 

0.00164659 -0.00104207 -0.00101912 

000635563 -0.000422485 -0.000398759 .. 

000251419 -0.000172685 -0.000159168 



0 .00164132 
,000662836 0 
.000269842 0 



0 

0. 

0.000110709 0.000101113 

]; 

end 

1 = length (g) ; 
f « 
f (1: 
f (1: 



zeros (1,2*1-1) ; 
2*1-1) = g; 

1-1) = reverse (g (2 : 1) ) ; 



end 



55 f = f ./ norm(f) ; 
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5 % 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 
10 % 



function [maxi , tree , T] = NonPrior2dStatTree (img, D, qmf , I , mult_sigma) 

15 [n,J] = quadlength(img) ; 
boxlen = n; 
boxcnt = 1 ; 
coef = img; 
coefl = I; 

20 

for d = 0 :D, 

for bx=0: (boxcnt -1) , 
for by=0: (boxcnt-1) , 

25 [lox, hix, loy, hiy] = quadbounds (d, bx, by , n) ; 

Quad = coef (lox: hix, loy : hiy) ; 
QuadI = coef I (lox : hix, loy :hiy) ; 

[si, s2] =size (Quad) ; 

30 

sigma=sqrt (1/ (sl*s2) * sum (sum (QuadI -* QuadI) ) ) ; 



% [y,estim] = NormNoise (Quad, qmf ) ; 

35 % sigma=l/estim 

T_node =mu 1 t_s i gma * s i gma / 

rl= abs (Quad) < T_node; 
40 r2=l-rl; 

Quad= sign (Quad) . * (abs (Quad) -sigma) ; 

pl=Quad. A 2; 

ptl»pl. *rl; 

p2=(T_node^2+sigma x 2) *ones (size (Quad) ) ; 
45 pt2=p2.*r2; 

tree(qnode (d,bx,by) ) = sum(sum(ptl) ) +sum (sum (pt2 ) ) ; 

T(qnode(d,bx,by) ) « T_node ; 

maxi (qnode (d,bx,by) ) =max (max (Quad) ) ; 

50 

% fin du bloc 

if d < D, % prepare for finer scales 
quad = DownQuad ( Quad , qmf , by , bx) / 
55 quadl o DownQuad (QuadI, qmf ,by,bx) ; 

coef (lox: hix, loy :hiy) = quad; 
coef I (lox: hix, loy :hiy) = quadl; 
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end 
end 
end 

boxlen = boxlen/2; 
boxcnt = boxcnt*2; 
end 



function [y,coef] = NormNoise (x, qmf ) 

% NormNoise Estimates noise level, Normalize signal to noise level 
10 1 

% Usage 

% [y,coef] = NormNoise (x, qmf ) 
% Inputs 

% x 1-d signal 

15 % qmf quadrature mirror filter 
% Outputs 

% y 1-d signal, scaled so wavelet coefficients 

% at finest level have median absolute deviation 1 . 

% coef estimation of 1/sigma 
20 % 

% Description 

% This is required pre-processing to use any of the DeNoising 

% tools on naturally-occurring data. 

% 

25 % See Also 

% WaveShrink, CPDeNoise , WPDeNoise 
% 

u = DownDyadHi (x, qmf ) ; 
s = median (abs (u) ) ; 
30 if s -= 0 

y m .6745 .* x ./s; 
coef = .6745 /s; 

else 

y = x; 

35 coef = 1; 

end 

% 

% Copyright (c) 1993-5. Jonathan Buckheit, David Donoho and Iain 
40 Johnstone 
% 

% Modified by Maureen Clerc and Jerome Kalifa, 1997 

% clerc@cmapx.polytechnique.fr, kalifa@cmapx .poly technique . fr 

% 

45 

% 

% Part of WaveLab Version 8 02 
% This is Copyrighted Material 
50 % For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 
% 



55 function ix = qnode (d, bx, by) 
% qnode -- Quad tree indexing 
% Usage 

-35- 



WO 01/84500 



PCT/US01/13737 



10 % 



20 



% ix = qnode (d,bx,by) 

% Inputs 

% d depth in splitting 

% bx,by box coordinates at that depth 

% Outputs 

% ix index of that node in tree 

% 

% Description 

% ix = 4*(d) + by + bx*2"(d) 



ix = 4 A (d) + by + bx*2^(d); 



15 % Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 



% 



function [lox, hix, loy , hiy] = quadbounds (d, bx, by, n) 
% quadbounds x, y bounds of quadlet 
% Usage 

25 % [lox, hix, loy, hiy] = quadbounds {d, bx, by , n) 

% Inputs 

% d depth in splitting 

% bx,by box coordinates at that depth 

% n extent of image 

30 % Outputs 

% lox lower x- index of quadlet 

% loy lower y- index of quadlet 

% hix upper x- index of quadlet 

% hiy upper y- index of quadlet 

35 % 

% Description 

% This routine is called by other, higher-level Wavelab functions 

% and is not intended to be useful for most users . 

% 

40 boxsize = n/ (2^d) ; 

lox = bx*boxsize+l; 
loy = by*boxsize+l; 
hix = (bx+1) *boxsize; 
hiy = (by+1) *boxsize; 

45 

% 

% Part of WaveLab Version 8 02 
% This is Copyrighted Material 
50 % For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 
% 



55 function [n,J] = quadlength (x) 

% quadlength Find length and dyadic length of square matrix 
% Usage 
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% [n,J] = quadlength (x) 
% Inputs 

% x 2-d image; size(n,n), n = 2^J (hopefully) 
% Outputs 
5 % n length (x) 

% J least power of two greater than n 
% 

% Side Effects 

% A warning message is issue if n is not a power of 2, 
10 % orifxis not a square matrix. 
% 

s = size (x) ; 
n m s (1) ; 
if s (2) ~= s (1) , 
15 disp ( 'Warning in quadlength: nr != nc') 

end 

k = 1 ; J = 0; while k < n , k=2*k; J = 1+J ; end ; 
if k -= n , 

disp ( 'Warning in quadlength: n != 2^ J') 

20 end 

% 

% Copyright (c) 1993. David L. Donoho 
% 

25 

% 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
30 % For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 
% 



35 functi on x = Sof tThresh (y , t) 

% Sof tThresh Apply Soft Threshold 
% Usage 

% x = Sof tThresh (y, t) 
% Input s 
40 % y Noisy Data 

% * t Threshold 
% Outputs 

% x sign(y) ( |y| -t)_+ 

% 

45 res = (abs (y) - t) ; 

res = (res -f abs (res)) /2; 
x = sign (y) . *res ; 

% 

50 % Copyright (c) 1993-5. Jonathan Buckheit, David Donoho and Iain 
Johnstone 
% 



55 % 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
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% For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 

% 

5 . 

function 

Inv_Recons=SPECTMacrTomoFiltback (num, preDenoise , postDenoise , mult_pre , 

mult_j?ost , f ileproj , typel , type2 ) 

% 

10 % Macro for tomographic reconstruction methods using SPECT Interfile 
data 

% With an adaptive selection of the wavelet packet basis 
% 

15 dime=128; 

sino = Decodef ile(f ileproj , num); 
sino ( : , : ) =Big2Little (sino (:,:)); 
Theta=0:2.82:3 60; % Projection angles 

20 

% Estimation of the Noise level 

qmf = MakeONFilter { ' Symmlet 1 , 4 ) ; 

25 [y,coefj = NormNoise (sino, qmf ) ; 
sigma=l/coef 

Inv_Degr=f iltback ( sino) ; 

30 if postDenoise==l 

% Uncomment if you want to use a Phantom Image as prior model for 

the 

% image 

35 % Phant=raw2mat ( 'phantoml2 8 .gray ' , 12 8, 12 8) / ^ 

B= (randn (size (sino) )) *sigma; % Noise generation 

Inv__B= filtback(B); 

40 D=3; 

% If one uses a Phantom to choose the basis 

% [maxi, tree,T] = Dis2dStatTree (Phant , D, qmf , Inv_B , mult_post) ; 

45 % if one uses the degraded data to choose the basis 

[maxi, tree, T] = 
NonPrior2dStatTree ( InvJDegr , D , qmf , Inv_B, mult_post) ; 

50 btree = Best^dBasis (tree, D) ; 

% ax = axis; hold; 

% Plot2dPartit ion (btree, ! y' ,ax,D) ; 

Rest_best=zeros (size (Inv_Degr) ) ; 
lim=2 ; 

55 

for i=l:lim, 
for j=l:lim, 
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[Inv__Degr trans] = cyclespin2 (Inv_Degr , i , j ) ; 

[Rest] =DiscrDenoise2Radon(btree, Inv__Degrtrans , 0*T,qmf , ' H 1 ) ; 
[Rest] = cyclespin2 (Rest, dime-i, dime- j ) ; 
Rest_best = Rest_best+Rest ; 
5 end 
end 

Rest_best=Rest_best/ (lim*lim) ; 
Inv_Recons=Rest_best ; 

10 end 



% 

15 % <c) Jerome Kalifa, 2000 
% kalifa@columbia.edu 
% 

function Res t_bes t=SPECTTomo (num, mult , f ileproj ,type) 

20 

% size of the images 
dime=128 ; 

25 % Read Interfile projections for the slice #num 
sino = Decodefile (f ileproj , num) ; 
% Big Endian Little Endian 

30 

sino( : , : ) =Big2Little (sino { :,:)); 
% Projection angles 
35 Theta=0 :2 . 82 :360; 

% Filter to be used by the wavelet (packet ) transform 
qmf « MakeONFilter ( 1 Symmlet 1 , 4) ; 

40 

% D= 2 or 3 is fine (3 slightly better, but slower) ; 
D=2; 

% Simulation of the backproj ected noise, 1st alternative 

45 

[y,coef] = NormNoise (sino, qmf ) ; 
sigma=l/coef ; 

B= (randn(size (sino) ) ) *sigma; 
50 % 2nd alternative, using numerical measurements 

%out=ThreshWave2 (sino, ' H 1 ) ; 
%B= sino -out ; 

55 % Generation of the backproj ected noise, and measurement of the 
amplitude of 
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% its coefficients in the WP tree. NB : If the noise B has been 
generated 

% with 1st alternative, it is possible to avoid this step by storing 
a array 

5 % A of noise amplification values for the WP coefficiens, and 
computing 

% AmpNoise=sigma. *A . 

Inv_B= filtback(B) ; 
10 AmpNoise= TomoStatTree (D, qmf , Inv_B) ; 

% Wavelet packet basis we want to use. 

if D==2 ; 
15 btree = zeros (1, 31) ; 

btree (1)=1; 
btree (5: 6) = [1 1] ; 
elseif D==3; 

btree = zeros (1 , 12 7 ) ; 
20 btree ( 1)=1 ; 

btree(5:7) = [l 11]; 
btree (26 :27)=[1 1] / 
btree (30: 31) = [1 1] ; 
end 

25 

% Reconstruction by BackProj ection, just ramp filtering 
Inv_Degr=f iltback (sino) ; 
30 % Initialisation of the final image 
Rest_best=zeros (dime, dime) ; 

% Translation-invariance. lim=4 slightly better than lim=2, but 
35 slower. 
lim=2 ; 

% Denoising in the wavelet packet tree. 

40 for i=l:lim, 

for j=l:lim, 

[Inv__Degr trans] = cyclespin2 (Inv_Degr , i , j ) ; 

[Rest] = TomoDenoise (btree, Inv__Degrtrans,AmpNoise, qmf, type, mult) ; 
[Rest] - cyclespin2 (Rest , dime-i , dime- j ) ; 
45 Rest_best = Rest_best+Rest ; 

end 
end 

Rest_best=Rest_best/ (lim*lim) ; 
50 Rest_Int=Rest_bes t ; 

% Denoising in the wavelet tree 

Rest_best=zeros (size (Rest_Int) ) ; 
55 btree=zeros (size (btree) ) ; 
btree (1)=1; 
btree (4) =1; 
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lim=2 ; 

for 1=1:11x11, 
for j=l :lim, 

5 [Rest_Inttrans] = cyclespin2 (Rest_Int, i, j ) ; 

[Rest] =TomoDenoise (btree , Rest_Inttrans , AmpNoise , qmf , 1 H 1 , mult) ; 

[Rest] = cyclespin2 (Rest,dime-i, dime-j) ,- 

Rest_best = Rest_best+Rest ; 
end 
0 end 



5 



0 



5 



0 



5 



0 



5 



0 



% Regularized reconstructed image 
Rest_best=Rest_best/ (lim*lim) ; 
% 

% (c) Jerome Kalifa, 2001 
% kalif aOcolumbia . edu 
% 



function out=ThreshWave2Af f (Noisy, type, TI , sigma, mult , L, qmf ) 

% ThreshWave2 Denoising of 2-d image with wavelet thresholding. 

% Usage 

% out=ThreshWave2 (Noisy, type,TI, sigma, mult , L, qmf) 
% Inputs 

% Noisy 2-d Noisy image, size (Noisy) = 2^J*2^J. 

% type for soft thresholding, , H I for hard 

thresholding . 

% Optional, Default=hard thresholding. 

% TI Enter 1 if you want translation-invariant 

denoising, 

0 if you don't. 

Optional , Def ault=non- invariant . 

Standard deviation of additive Gaussian White 



% 
% 

% sigma 
Noise . 
% 

filtering . 
% 

% mult 

threshold. 

% 

% 

% L 

% 
% 

% qmf 
% 

% Outputs 
% out 
% 
% 



Enter 0 if you want sigma to be estimated by median 

Optional, Def ault=estimated by median filtering. 

Multiplier of sigma to obtain the value of the 

Optional, Default = sqrt (2 *log (n) ) , where n is the 
number of pixels/ 

Low- Frequency cutoff for shrinkage (e.g. L«4) 
Should have L << J I 
Optional, Default = 3. 

Quadrature Mirror Filter for Wavelet Transform 
Optional, Default = Symmlet 4. 

Estimate, obtained by applying thresholding on 
wavelet coefficients. 



n=length (Noisy) ; 
if nargin < 7 , 

qmf = MakeONFilter ( 1 Symmlet ' ,4) ; 

end 

if nargin < 6, 
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Ii = 3; 

end 

if nargin < 5/ 

mult = sqrt (4*log(n) ) ; 

5 end 

if nargin < 4 , 

[y,coef] = NormNoise (Noisy, qmf) ; 
sigma=l/ coef ; 

end 

10 if nargin < 3, 

TI = 0; 

end 

if nargin < 2 , 
type = »H' ; 

15 end 

if sigma==0, 

[y,coef] = NormNoise (Noisy, qmf ) ; 
sigma=l/ coef ; 

20 end 

thresh= sigma*mult ; 
out=zeros (l,n) ; 



if TI==1, 

if strcmp (type, 1 S 1 ) / 

out=TIDenoiseSof t2 (Noisy , L, qmf , thresh) ; 
end 

30 if strcmp (type, 'H 1 ), 

out=TIDenoiseHard2 (Noisy , L, qmf , thresh) ; 
end 

if strcmp (type, ' A* ) , 

out=TIDenoiseAf f ine2 (Noisy, L, qmf , thresh) ; 
35 end 
else 

if strcmp (type, 'S ') , 

out-ST2 (Noisy, L, qmf , thresh) ; 

end 

40 if strcmp (type, 'H» ) , 

out=HT2 (Noisy, L, qmf, thresh) / 

end 

if strcmp (type, 'A 1 ) , 

out=AT2 (Noisy, L, qmf , thresh) ; 

45 end 
end 

% Written by Maureen Clerc and Jerome Kalifa, 1997 

% clerc@cmapx . polytechnique . f r , kalif a@cmapx . polytechnique . f r 

50 



55 % 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
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% For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 

% 



TomoDenoise (basis , img, AmpNoise , qmf , type , mult ) 



function out 
% 

% Usage 

% out = TomoDenoise (basis, img, T, qmf ,D,type) 
10 % 

[n,J] = quadlength(img) ; 
coef = img; 

% 

15 % initialize tree traversal stack 
% 

stack « zeros (3,4* (J+l) +1) ; 
k = 1; 

stack(:,k) = [0 0 0]'; % d, bx, by 



20 % 



25 



30 



35 



40 



45 



55 



while (k > 0) , 



thresh) ; 



50 tnresh) ; 



thresh) ; 



% pop stack 
d = stack (l,k) ; 
bx « stack (2, k) ; 
by = stack (3 , k) ; 
k=k-l ; 

[lox hix loy hiy] « quadbounds (d, bx, by, n) / 
Quad = coef (lox: hix, loy : hiy) ; 

if (basis (qnode(d,bx, by)) — 0) , % nonterminal node 
quad = DownQuad (Quad , qmf , bx , by) ; 
coef (lox:hix, loy :hiy) = quad; 

% push stack 

k = k+1; stack(:,k) - [(d+1) (2*bx) (2*by) ] 

k = k+1; stack(:,k) = [(d+1) (2*bx+l) (2*by) ] 

k - k+1; stack(:,k) = [(d+1) (2*bx) (2*by+l) ] 

k = k+1; stack (:,k) = [(d+1) (2*bx+l) (2*by+l) ] 



else 



coef f_max=max (max (abs (Quad) ) ) ; 
if (coeff_max > (6 . *AmpNoise) ) , 

thresh=mult . *AmpNoise (qnode (d,bx,by) ) ; 

if strcmp(type, f A f ) , 

coef (lox: hix, loy: hiy) = Af f ineThresh (Quad, 

end 

if strcmp(type, 'S 1 ) , 

coef (lox: hix, loy: hiy) = Sof tThresh (Quad, 

end 

if strcmp(type, 'H f ) , 

coef (lox: hix, loy: hiy) = HardThresh (Quad, 

end 
else 

% coeff max 
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coef (lox: hix, loy : hiy) = zeros (hix-lox+1 , hiy- 

loy+1) ; 

end 

end 

end 

out = IPT2_WP (basis, coef , qmf ) ; 
function AmpNoise = TomoS tat Tree (D, qmf , I) 



10 



15 



35 



[n,J] = quadlength (I) ; 
boxlen = n; 
boxcnt = 1 ; 
coefl = I; 

for d = 0:D, 

for bx=0 : (boxcnt-1) , 
for by=0: (boxcnt-1) , 

20 [lox, hix, loy, hiy] = quadbounds (d, bx, by , n) ; 

QuadI = coef I (lox: hix, loy : hiy) ; 

[si, s2] =size (QuadI) ; 

AmpNoise_node= sqrt (1/ (sl*s2 ) *sum (sum (Quadl . *QuadI) ) ) ; 
25 AmpNoise (qnode (d,bx,by) ) = AmpNoise_node ; 

% prepare for finer scales 
if d < D, 

30 quadl = DownQuad (QuadI, qmf ,by,bx) ; 

coef I (lox: hix, loy : hiy) = quadl; 
end 
end 
end 



boxlen = boxlen/" 2 ; 
boxcnt = boxcnt*2; 
end 



40, function Q = UpQuad (quad, qmf , xbit , ybit) 

% UpQuad Merge four subbands into 2d image 

% Usage 

% Q = UpQuad (quad, qmf , xbit, ybit j 

% Inputs 

45 % Quad quadlet to be joined 4->l channel 

% qmf orthogonal quadrature mirror filter 

% xbit x-position of quadlet in 2-d freq plane 

% ybit y-position of quadlet in 2-d freq plane 

% Outputs 

50 % q quadlet after splitting 
% 

% Description 

% A 2-d signal split into four channels: (Lo_x,Lo_y) , 

% (IiO_x,Hi_y) , (Hi_x,Lo_y) , (Hi__x,Hi_y) is unsplit 

55 % into original image. 
% 

% This routine is called by other, higher-level WaveLab 
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10 % 



25 



% functions and may not be useful for many users. 
% 

% See Also 

% IPT2_WP, DownQuad 
% 

[nr,nc] = size(quad); 
Q = zeros (nr,nc) ; 
ybit = rem(ybit,2); 
xbit = rem(xbit , 2) ; 

bot = 1 : (nr/2) ; 
top = bot + (nr/2) ; 



for j=l:nc, 
15 if ybit, 

lpchan = quad (top , j ) 1 ; 
hpchan - quad (bot , j ) ' ; 

else 

lpchan = quad (bot , j ) 1 ; 
20 hpchan = quad (top, j) f ; 

end 



Q ( : , j ) = UpDyadLo (lpchan, qmf ) ■ + UpDyadHi (hpchan, qmf) ' 

end 

bot = 1: (nc/2) ; 
top = bot + (nc/2) ; 



for i = 1 : nr , 
30 if xbit, 

lpchan ■ Q(i,top); 

hpchan = Q(i,bot); 

else 

lpchan = Q(i,bot); 

35 hpchan = Q(i,top); 

end 



40 



Q (i , : ) = UpDyadLo (lpchan, qmf) + UpDyadHi (hpchan, qmf) ; 

end 



% 

% Part of WaveLab Version 802 
45 % This is Copyrighted Material 

% For Copying permissions see COPYING . m 

% Comments? e-mail wavelab®stat . Stanford . edu 
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APPENDIX B 

function dwt = FWTJVTrou (x, L) ; 
% 

5 %' FWT_ATrou -- Fast Dyadic Wavelet Transform (periodized, orthogonal) 
% Usage 

% dwt = FWT_ATrou(x,L) 
% Inputs 

% x 1-d signal; length (x) = 2"J = n 

10 % L Coarsest Level of V_0 ; L << J 

% Outputs 

% dwt an n times J-L+l matrix 

% giving the wavelet transform of x at all dyadic scales. 

% 

15 % Description 

% To reconstruct use IWT_ATrou 
% 

% See Also 

% IWT_ATrou, Make ATrouFi Iter 
20 % 

[lodyadf , dlodyadf , hidyadf , dhidyadf ] = MakeATrouFilter ( 1 Spline ' , 3 ) ; 

[n,J] - dyadlength (x) ; 
25 D = J-L; 

dwt = zeros (n,D+l) ; 
x = Shape As Row (x) ; 
dwt (1,1) = x 1 ; 

30 for d = 1:D, 

s = dwt ( : , 1) 1 ; 
S2 = S; 

for j = 1:2" (d-l) 
s2 = lshif t (s2) / 
35 end 

dwt(:,d+l) = iconv (hidyadf , s) ' ; 
for j = 1:2" (d) 

p = lshif t (dwt (: ,d+l) ') ; 
dwt ( : , d+1) = p' ; 
40 end 

dwt(:,l) = iconv (lodyadf , s2) ' ; 

f = zeros (1, 2 * length (lodyadf) ) ; 

f (1:2 :2*length( lodyadf ) -1) = lodyadf; 

45 

f2 = zeros (1, 2*length (hidyadf )) ; 

f 2 (1:2 :2*length (hidyadf) -1) = hidyadf ; 

lodyadf = f; 
50 hidyadf = f2; 

end; 

% Written by Maureen Clerc and Jerome Kalifa, 1997 

% clerc@cmapx.polytechnique.fr, kalifa@cmapx.polytechnique.fr 

55 
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% 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
5 % Comments? e-mail wavelab@stat.stanford.edu 
% 



function s = I WT_ATr ou ( dwt , L) ? 
% 

10 % IWT_ATrou -- Inverse Dyadic Wavelet Transform 
% Usage 

% s = IWT_ATrou (dwt,L) 
% Inputs 

% dwt an n times J-L+l matrix 
15 % L Coarsest Level of V_0; L << J 
% Outputs 

% s original 1-d signal; length (x) = 2" J = n 

% Description 

% 1. filters are obtained with MakeATrouFilter 
20 % 2. usually, length (qmf) < 2" (L+l) 

% 3. The transformed signal can be obtained by FWT_ATrou 
% See Also 

% FWT_ATrou , MakeATr ouFi Iter 
% 

25 [n,b] = size (dwt ); 
J=b+L-1; 
if - (2" J==n) 

fprintf ( 1 problem in matrix dimensions'); 
end 

30 s = dwt ( : ,1) • ; 
D = b-1; 

[lodyadf , dlodyadf , hidyadf , dhidyadf ] « MakeATrouFilter ( ' Spline 1 , 3 ) ; 

f = zeros (1,2" (D-l) *length (dlodyadf )) ; 
35 f (1:2" (D-l) :2" (D-l) *length (dlodyadf) +1-2" (D-D) = dlodyadf; 

f2 = zeros (1,2" (D-l) *length (dhidyadf )) ; 

f2 (1:2" (D-l) : 2" (D-l) *length (dhidyadf ) +1-2" (D-l) ) = dhidyadf; 

40 for d= D-l: -1:0 

for j « 1:2 . *2" (d) , 2 

s = Ishift(s) ; 
end 

for j m 1:2" (d+1) 
45 p = rshift (dwt ( : ,d+2) ' ) ; 

dwt ( : , d+2 ) = p ' ; 
end 

for j =1:3 . *2" (d) , 
50 p=lshif t (dwt { : , d+2) ' ) ; 

dwt ( : , d+2 ) =p 1 ; 
end 

s=0.5* (iconv(f , s) +iconv(f 2,dwt ( : ,d+2) ■ ) ) ; 

55 f .= zeros (1,2" (d-l) *length ( dlodyadf) ) ; 

f (1:2" (d-l) : 2" (d-l) *length (dlodyadf) +1-2" (d-l)) = dlodyadf; 
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f2 = zeros (1,2^ (d-1) * length (dhidyadf )) ; 

f2 (1:2" (d-1) :2 A (d-1) * length (dhidyadf ) +1-2" (d-1) ) = dhidyadf; 
end 

5 % Written by Maureen Clerc and Jerome Kalifa, 1997 

% clerc@cmapx.polytechnique . f r, kalif a@cmapx.polytechnique . f r 



% 

10 % Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 
% 

15 



function y = lshif t (x) 

% lshif t Circular left shift of 1-d signal 
% Usage 
20 % 1 = lshif t(x) 
% Inputs 

% x 1-d signal 
% Outputs 
% 1 1-d signal 
25 % l(i) = x(i+l) except 1 (n) = x(l) 

% 

y = [ x( 2:length(x) ) x(l) ]; 

30 % 

% Copyright (c) 1993. Iain M. Johnstone 
% 



35 % 

% Part of WaveLab Version 802 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 
40 % 



function [lodyadf , dlodyadf , hidyadf /dhidyadf]. = 
MakeATrouFilter (Type , Par) 
45 % 

% MakeATrouFilter Generate Biorthonormal Quadratic Spline Filter 

Pair 

% 

% Usage 

50 % [lodyadf ,dlodyadf,hidyadf, dhidyadf] = MakeATrouFilter (Type , Par) 
% Inputs 

% Type string, one of: 
% 'Spline' 
% Par Par = 3 only 
55 % 

% Outputs 

% lodyadf low-pass dyadic filter 
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% dlodyadf dual low-pass dyadic filter . 
% hidyadf high-pass dyadic filter 

% dhidyadf dual high-pass dyadic filter 
sqrt2 = sqrt (2) ; 
5 if strcmp (Type, 'Spline' ) , 

lodyadf = [0.125 0.375 0.375 0.125] . *sqrt2; 

% for Haar 

% lodyadf « [0. 0.5 0.5 0].*sqrt2; 

dlodyadf = [0.125 0.375 0.375 0 . 125] . *sqrt2 ; 
10 hidyadf = [.5 -.5] . *sqrt2; 

dhidyadf = [-.03125 -.21875 -.6875 .6875 .21875 
.03125] . *sqrt2; 
end 

15 % Written by Maureen Clerc and Jerome Kalifa, 1997 

% clerc@cmapx.polytechnique . f r , kalif a@cmapx .poly technique . f r 



% 

20 % Part of WaveLab Version 8 02 
% This is Copyrighted Material 
% For Copying permissions see COPYING. m 
% Comments? e-mail wavelab@stat.stanford.edu 
% 

25 



function row = ShapeAsRow(sig) 
% ShapeAsRow Make signal a row vector 
% Usage 
30 % row = ShapeAsRow (sig) 
% Inputs 

% sig a row or column vector 

% Outputs 

% row a row vector 

35 % 

% See Also 

% ShapeLike 

% 

row = sig ( : ) 1 ; 



% 

% Part of WaveLab Version 8 02 
% This is Copyrighted Material 
45 % For Copying permissions see COPYING. m 

% Comments? e-mail wavelab@stat.stanford.edu 
% 
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CLAIMS 

1 . A method for reconstructing image data corresponding to a desired 
two-dimensional image from initial tomographic data associated with a plurality of 
initial tomographic projections of said desired image, comprising: 

backprojecting said initial tomographic data to form backprojected 
image data; 

selecting a vector family to decompose said backprojected image data 
from the group consisting of 1-D orthogonal wavelet packet basis, 2-D 
orthogonal wavelet packet basis, 1-D biorthogonal wavelet packet basis, 2-D 
biorthogonal wavelet packet basis, 1-D wavelet packet frame, 2-D wavelet 
packet frame, 1-D complex wavelet frame, 2-D complex wavelet frame, 1-D 
complex wavelet packet frame, and 2-D complex wavelet packet frame; 

decomposing said backprojected image data to generate decomposed 
coefficients associated with said selected vector family; 

diagonally operating on said decomposed coefficients to generate 
modified decomposed coefficients; and 

recomposing said image data corresponding to said desired two- 
dimensional image from said modified decomposed coefficients. 

2. A method for reconstructing image data corresponding to a desired 
20 two-dimensional image from initial tomographic data associated with a plurality of 

initial tomographic projections of said desired image, comprising: 

selecting a vector family to decompose said initial tomographic data 
from the group consisting of 1-D orthogonal wavelet packet basis, 2-D 
orthogonal wavelet packet basis, 1-D biorthogonal wavelet packet basis, 2-D 
25 biorthogonal wavelet packet basis, 1-D wavelet packet frame, 2-D wavelet 

packet frame, 1-D complex wavelet frame, 2-D complex wavelet frame, 1-D 
complex wavelet packet frame, and 2-D complex wavelet packet frame; 

decomposing said initial tomographic data into decomposed 
coefficients associated with said selected vector family; 
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diagonally operating on said decomposed coefficients to generate 
modified decomposed coefficients; and 

recomposing said modified decomposed coefficients to generate 
regularized data associated with a plurality of regularized tomographic 
5 projections; and 

backprojecting said regularized data to form said image data 
corresponding to said desired two-dimensional image. 

3. A method for reconstructing image data corresponding to a desired 
two-dimensional image from initial tomographic data associated with a plurality of 
initial tomographic projections of said desired image, comprising: 

selecting a vector family to decompose said initial tomographic data 
from the group consisting of 1-D orthogonal wavelet packet basis, 2-D 
orthogonal wavelet packet basis, 1-D biorthogonal wavelet packet basis, 2-D 
biorthogonal wavelet packet basis, 1-D wavelet packet frame, 2-D wavelet 
packet frame, 1-D complex wavelet packet frame, and 2-D complex wavelet 
packet frame; 

decomposing said initial tomographic data into decomposed 
coefficients associated with said selected vector family; 

diagonally operating on said decomposed coefficients to generate 
modified decomposed coefficients; 

selectively amplifying said modified decomposed coefficients to 
generate selectively amplified coefficients; 

recomposing said selectively amplified coefficients to generate 
regularized data associated with a plurality of regularized tomographic 
projections; and 

radially interpolating said regularized data to form said image data 
corresponding to said desired two-dimensional image. 

4. A method for reconstructing image data corresponding to a desired 
two-dimensional image from initial tomographic data associated with a plurality of 

30 initial tomographic projections of said desired image, comprising: 
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selecting a vector family to decompose said initial tomographic data 
from the group consisting of 1-D orthogonal wavelet packet basis, 2-D 
orthogonal wavelet packet basis, 1-D biorthogonal wavelet packet basis, 2-D 
biorthogonal wavelet packet basis, 1-D wavelet packet frame, 2-D wavelet 
5 packet frame, 1-D complex wavelet packet frame, and 2-D complex wavelet 

packet frame; 

decomposing said initial tomographic data into decomposed 
coefficients associated with said selected vector family; 

selectively amplifying said decomposed coefficients to generate 
10 selectively amplified coefficients; 

diagonally operating on said selectively amplified coefficients to 
generate modified decomposed coefficients; 

recomposing said modified decomposed coefficients to generate 
regularized data associated with a plurality of regularized tomographic 
15 projections; and 

radially interpolating said regularized data to form said image data 
corresponding to said desired two-dimensional image. 

5. The method of claim 1, 2, 3 or 4 wherein said diagonally operating 
comprises thresholding the data operated on. 

20 6. The method of claim 1, 2, 3 or 4 wherein said selecting a vector family 

step further comprises: 

deriving a model of noise present in said data to be decomposed; 
repeatedly selecting a wavelet packet decomposition from a 
predetermined plurality of wavelet packet decompositions; 
25 calculating an estimation of final estimation error based on at least said 

selected wavelet packet decomposition and said noise model; and 

choosing a wavelet packet decomposition that minimizes the 
estimation of final estimation error calculated in said repeatedly selecting step. 
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7. The method of claim 6 wherein said calculating an estimation of final 
estimation error is also based on prior information about said desired two-dimensional 
image. 

8. A method for reconstructing image data corresponding to a plurality of 
5 two-dimensional slices of a desired three-dimensional image from initial tomographic 

data associated with a plurality of tomographic projections of said plurality of two- 
dimensional slices, comprising: 

selecting a vector family to decompose said initial tomographic data 
selected from the group consisting of wavelet vector family and wavelet 
1 0 packet vector family; 

decomposing said initial tomographic data to generate decomposed 
coefficients associated with said selected vector family; 

diagonally operating on said decomposed coefficients to generate 
modified decomposed coefficients; and 
1 5 recomposing said modified decomposed coefficients to generate 

regularized data associated with a plurality of regularized tomographic 
projections of said plurality of two-dimensional slices of said desired image; 
and 

backprojecting said regularized data to form said image data 
20 corresponding to said plurality of two-dimensional slices of said desired three- 

dimensional image. 

9. A method for reconstructing image data corresponding to a plurality of 
two-dimensional slices of a desired three-dimensional image from initial tomographic 
data associated with a plurality of tomographic projections of said plurality of two- 

25 dimensional slices, comprising: 

backprojecting said initial tomographic data to form backprojected 
image data associated with said plurality of two-dimensional slices; 

selecting a wavelet vector family to decompose said backprojected 
image data; 
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decomposing said backprojected image data to generate decomposed 
coefficients associated with said selected vector family; 

diagonally operating on said decomposed coefficients to generate 
modified decomposed coefficients; and 
5 recomposing said modified decomposed coefficients to form said 

image data corresponding to said plurality of two-dimensional slices of said 
desired three-dimensional image. 

10. The method of claim 8 further comprising performing slice-by-slice 
two-dimensional regularization on said initial tomographic data to form regularized 

10 data corresponding to a series of regularized two-dimensional slices before said 

decomposing step wherein said decomposing step is performed on the entirety of said 
regularized data. 

1 1 . The method of claim 9 further comprising performing slice-by-slice 
two-dimensional regularization on said initial tomographic data to form regularized 

15 data corresponding to a series of regularized two-dimensional slices before said 

backprojecting step wherein said backprojecting step is performed slice-by- slice on 
the entirety of said regularized data 

12. The method of claim 9 further comprising performing slice-by-slice 
two-dimensional regularization on said backprojected image data to form regularized 

20 data corresponding to a series of regularized two-dimensional slices before said 

decomposing step wherein said decomposing step is performed on the entirety of said 
regularized data 

13 . The method of claim 8 or 9 wherein said diagonally operating 
comprises thresholding the data operated on. 

25 14. The method of claim 8 or 9 wherein said selecting a vector family step 

further comprises: 

deriving a model of noise present in said data to be decomposed; 
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repeatedly selecting a wavelet packet decomposition from a 
predetermined plurality of wavelet packet decompositions; 

calculating an estimation of final estimation error based on at least said 
selected wavelet packet decomposition and said noise model; and 
5 choosing a wavelet packet decomposition that minimizes the 

estimation of final estimation error calculated in said repeatedly selecting step. 

15. The method of claim 14, wherein said calculating an estimation of final 
estimation error is also based on prior information about said plurality of two 
dimensional slices of said desired three-dimensional image. 
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